Nonsequential two-photon double ionization of helium 



J. Feist, 1 '* S. Nagele, 1 R. Pazourek, 1 E. Persson, 1 B. I. Schneider, 2, 3 L. A. Collins, 4 and J. Burgdorfer 1 

1 Institute for Theoretical Physics, Vienna University of Technology, 1040 Vienna, Austria, EU 
2 Physics Division, National Science Foundation, Arlington, Virginia 22230, USA 
3 Electron and Atomic Physics Division, National Institute of Standards and Technology, Gaithersburg, Maryland 20899, USA 
J * Theoretical Division, T-4, Los Alamos National Laboratory, Los Alamos, New Mexico 87545, USA 

(Dated: May 7, 2008) 

We present accurate time-dependent ab initio calculations on fully differential and total inte- 
grated (generalized) cross sections for the nonsequential two-photon double ionization of helium 
at photon energies from 40 to 54 eV . Our computational method is based on the solution of the 
time-dependent Schrodinger equation and subsequent projection of the wave function onto Coulomb 
waves. We compare our results with other recent calculations and discuss the emerging similarities 
and differences. We investigate the role of electronic correlation in the representation of the two- 
electron continuum states, which are used to extract the ionization yields from the fully correlated 
final wave function. In addition, we study the influence of the pulse length and shape on the cross 
sections in time-dependent calculations and address convergence issues. 

PACS numbers: 32.80.Rm, 32.80.Fb, 42. 50. Hz 



INTRODUCTION 



Double ionization of helium has long been of great in- 
terest in atomic physics since it provides fundamental 
insights into the role of electronic correlation in the full 
three-body Coulomb break-up process. Understanding 
the dynamics in this simple, two-electron system is cru- 
cial to understanding more complex atoms and even sim- 
ple molecules [1-16]. With the advent of intense light 
sources in the vuv and xuv region [17-23], the focus of 
interest has switched from one-photon to multi-photon 
processes. Specifically, two-photon double ionization has 
been the subject of intense theoretical studies. Several 
authors [24-33] have calculated generalized cross sections 
for the nonsequential two-photon ionization process in 
the energy range from 40 to 54 eV using a wide variety 
of computational methods. Despite considerable efforts, 
quantitative agreement between the different calculations 
has not yet been reached. The reasons for the remain- 
ing discrepancies are the subject of ongoing discussions. 
In particular, there have been speculations that the rep- 
resentation of the double continuum might be responsi- 
ble for the existing differences. Nevertheless, even for 
methods which take correlation into account in the fi- 
nal states, the cross sections obtained still disagree, and 
a systematic change in the results due to the improved 
treatment of electronic correlation has not yet been ob- 
served. Recently, experimental data has become avail- 
able as well. Hasegawa, Nabekawa et al. [34, 35] used the 
27th harmonic at 41.8 eV of a femtosecond pulse from a 
Thsapphire laser, and Sorokin et al. [36] performed their 
experiment at the FLASH free-electron laser in Hamburg 
at 42.8 eV. While this is a good beginning, the uncer- 
tainties in the data are too large to help in resolving 



differences in the theoretical calculations. With the cur- 
rent rapid progress in intense xuv sources, further exper- 
iments that cover larger energy ranges can be expected 
in the future which might help to clarify the situation. 

Calculations for two-photon ionization employ either 
a time-independent (TI) or a time-dependent (TD) ap- 
proach. TI methods involve either lowest-order pertur- 
bation theory (LOPT) or i?-matrix Floquet theory. TD 
methods are based on a direct solution of the time- 
dependent Schrodinger equation and are therefore not 
restricted to any given order of the perturbation. Thus, 
they can be applied equally well to the strong field 
regime. In the present case of moderate intensities of the 
xuv field (~ 10 12 W/cm 2 ), corrections to LOPT are ex- 
pected to be small. The decisive advantage of TD meth- 
ods comes here from a different aspect. Namely, TI cal- 
culations of processes involving correlated two-electron 
final states in the continuum, 4 , ki,k 2 ( r i> r 2), require the 
knowledge of the final state in the entire configuration 
space in order to calculate the transition amplitude 
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where C/W is the transition operator for an TV-photon 
process (N = 2 in the following). As the numerical or 
analytical determination of accurate correlated contin- 
uum final states remains a challenge, evaluation of Eq. 1 
involves, inevitably, additional approximations that are 
difficult to control. Adding the time as an additional de- 
gree of freedom to the six spatial dimensions of the two- 
electron problem allows one to bypass the determination 
of ^ki,k 2 - Propagating the wave packet for sufficiently 
long times enables the extraction of the relevant dynami- 
cal information entirely from the asymptotic region where 
electron correlations become negligible. Moreover, resid- 
ual errors can be controlled by systematically varying 
the propagation time. This advantage comes along with 
a distinct disadvantage: Results will, in general, depend 
on the time-structure imposed on the external perturba- 
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tion, specifically on the duration and temporal shape of 
the xuv pulse. A comparison with TI calculations on the 
level of (generalized) cross sections therefore requires a 
careful extraction of information and checks of the inde- 
pendence from pulse parameters. 

In our theoretical approach, we solve the time- 
dependent Schrodinger equation (TDSE) using the time- 
dependent close-coupling (TDCC) scheme, cf. [24-26, 
37]. For the spatial discretization, we employ a finite 
element discrete variable representation (FEDVR), and 
the temporal propagation is performed by means of the 
short iterative Lanczos (SIL) procedure with adaptive 
time-step control. We present detailed convergence tests 
as a function of pulse duration, pulse shape, duration 
of propagation, gauge, spatial grid structure, and partial 
wave decomposition. For two-photon double ionization 
in the photon energy range 40 to 50 eV, we reach an ac- 
curacy on the 2% level. Our present results are compared 
with available experimental and theoretical data. Atomic 
units are used unless indicated otherwise. 



II. METHOD OF PROPAGATION 

The interaction of a helium atom (with infinite nuclear 
mass) with linearly polarized light is described by the 
Hamiltonian 
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where the interaction with the electromagnetic field in 
the dipole approximation is either given in length gauge 
by 

Ht m = E(t) (z 1 + z 2 ) (3) 
or in velocity gauge by 

— (p,,i + P,,2 + —j 1 • 4 
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If the exact solution were available, the two gauges would 
be strictly equivalent. Within approximate solutions, 
however, discrepancies may arise. The degree of gauge 
dependence can therefore be exploited as a measure for 
the convergence of the numerical solution toward the ex- 
act solution. 



A. Time-dependent close-coupling 

In order to solve the time-dependent Schrodinger equa- 
tion 
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we expand the six-dimensional wave function ^(ri, r 2 ) in 
coupled spherical harmonics 
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Substitution of Eq. 6 into Eq. 5 yields a system of cou- 
pled partial differential equations in (r\, r 2 ,t) , the TDCC 
equations [37] 
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where in practice the sums have to be truncated at cer- 
tain maximum angular momenta (L m ax, h,max, ^2,max) ■ 

We restrict ourselves to partial waves with total M = 
since M is conserved for linearly polarized laser fields. 
The helium singlet (S = 0) ground state is space- 
symmetric. As the spin quantum number S is con- 
served in the dipole approximation, the wave func- 
tion is space-symmetric for all times, implying that 
Rf 2 i 1 (i" 2 ,ri,t) = Rf i (ri,r 2 ,t) . As a consequence, the 
functions Rf l2 (ri,r 2 ,t) need not be stored for l 2 > h- 
Together with the parity selection rules, which only allow 
certain combinations (L,l\,l 2 ), this greatly reduces the 
numerical effort for solving the close-coupling equations. 



B. Spatial discretization 

For the discretization of the radial functions 
Rf l2 (ri,r 2 ,t), we employ a finite element discrete vari- 
able representation (FEDVR) [38-41]. We divide the ra- 
dial coordinates into finite elements in each of which the 
functions Rf l2 are represented in a local DVR basis with 
a corresponding Gauss-Lobatto quadrature to ensure the 
continuity of the wave function at the element bound- 
aries. This method leads to sparse matrix representations 
of the differential operators and to a diagonal potential 
matrix (within quadrature accuracy). Additionally, the 
boundary condition at r± = r 2 = can be easily fulfilled 
by omitting the first basis function in the first finite el- 
ement. The derivative discontinuity in the partial wave 
expansion of the electron-electron interaction at r± = r 2 , 
on the other hand, demands special treatment to guar- 
antee an accurate representation of the Hamiltonian in 
the FEDVR basis [42, 43]. 



C. Temporal propagation 

For the temporal propagation of the solution of the 
coupled equations (8), we employ the short iterative 
Lanczos (SIL) method [44-46] with adaptive time-step 
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control. The initial He ground state is obtained by re- 
laxing an arbitrary test function in imaginary time. In 
the SIL method, sometimes also referred to as Arnoldi- 
Lanczos algorithm, the time evolution operator 



U(t, t + At) ~ exp -iH(t)At 
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is represented by an N x N matrix UW in an N- 
dimensional Krylov subspace which is formally generated 
by repeated action of H on an initial state \&(t) and sub- 
sequent Gram-Schmidt orthogonalization. 

The Lanczos algorithm is very effective because the 
matrix H^ is tridiagonal and can be directly obtained 
by use of a three-term recursion relation involving mainly 
matrix- vector and scalar products. The sparse structure 
of the kinetic energy matrices due to the division of co- 
ordinate space into finite elements enables efficient par- 
allelization. In our calculations, which have primarily 
employed computers based on cluster architectures, we 
have observed linear scaling of the computational speed 
up to 450 processors. This gives us the possibility to 
employ pulses with comparably long durations (up to a 
few femtoseconds) in our simulations. In addition, exten- 
sive numerical convergence tests can be performed within 
reasonable time (cf. section HID). 

The SIL method also allows for a convenient error 
control since for a given subspace order N the dif- 
ference of the propagated wave function *W(t + At) 
to the lower-order approximation \I>( iv_1 )(f + At) can 
be calculated with only little extra effort and may 
be used as a tolerance parameter during propagation. 
In our implementation, we use a Krylov subspace of 
fixed order (with TV = 12 for all calculations pre- 
sented in this work) and an adaptive time step so 
that ||#W(f + At) - ^ (N -V{t + At) || 2 is smaller than 
a given tolerance parameter (typically in the order of 
10~ 20 ). Hence, we achieve a high-order approximation of 
the exponential function (Eq. 9) , and the temporal prop- 
agation is explicitly unitary and unconditionally stable. 
In addition, high accuracy is guaranteed due to the au- 
tomatic adjustment of the time step. 



III. EXTRACTING DYNAMICAL 
INFORMATION 

Subsequent to the time propagation, the information 
on excitation and ionization probabilities and on differ- 
ential and total cross sections must be extracted from 
the wave packet that is represented on a grid in a finite 
domain in coordinate space. This is a non-trivial task as 
straightforward projection is, for all practical purposes, 
precluded. Ideally, one would project onto asymptotic 
eigenstates of He, including its single and double ioniza- 
tion continua. As no closed solutions for the double con- 
tinuum are available, this is not directly possible. Numer- 
ical diagonalization of the Hamiltonian is not feasible ei- 
ther. First, the basis is too large to allow diagonalization. 



Second, the inclusion of the correct asymptotic boundary 
conditions in a numerical solution is highly non-trivial 
and computationally expensive. In addition, avoiding 
the construction of these eigenstates is one of the sig- 
nificant advantages of the TD approach. We circumvent 
the need for projection onto exact eigenstates by exploit- 
ing the fact that we can propagate the wave packet for 
long times after the conclusion of the pulse. Thereafter, 
the three-body Coulomb system (nucleus and two elec- 
trons) has reached (near asymptotically) large distances 
such that projection onto approximate energy eigenstates 
is possible without significant error. 

The states we eventually use to represent this contin- 
uum are eigenstates of the He Hamiltonian without the 
electron-electron interaction H12 = |fi — ?2 , i.e., prod- 
ucts of two Coulomb waves with Z = 2. This amounts 
to assuming that the electrons are far enough apart that 
their influence on each other can be neglected. By the 
same token, products of plane waves would be equally 
applicable when the electrons are far from the nucleus 
and its Coulomb potential could therefore be neglected. 
The advantage of using Coulomb waves is that orthogo- 
nality to bound states is built in. As the singly ionized 
part of the wave function can to a very good approxima- 
tion be written as the product of an unperturbed bound 
state of He + and a Coulomb wave (with effective charge 
Z = 1), the projection onto products of Coulomb waves is 
automatically orthogonal to the singly ionized part and 
no additional screening of the wave function has to be 
performed. Neglecting the electron-electron interaction, 
which is purely repulsive, introduces a small energy shift 
in the spectrum that decreases as the distance between 
the electrons increases. This effect can be controlled by 
varying the time of projection. 

We obtain the energy-space wave function by project- 
ing the spatial wave function onto products of energy- 
normalized Coulomb waves 4>e,i, i.e., 
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Switching from the angular momentum representa- 
tion {L,l\,l 2 ) to the representation in angular vari- 
ables (fli, fia), the six-dimensional (effectively, five- 
dimensional) distribution can be written as [26, 47] 
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where ai = argr(l+/+i77) are the Coulomb phases. From 
Eq. 11, reduced probability distributions can be deter- 
mined by integrating over unobserved degrees of freedom. 
For example, integrating over the solid angles (f2i,f2 2 ) 
gives the energy distribution (E 1 ,E 2 ) of the electron pair 
(Fig. 1). 
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FIG. 1: Energy distribution after two-photon double ioniza- 
tion from two different laser pulses with a mean energy of 
(uj) — 42 eV. The pulse (a) has a sin 2 envelope of total du- 
ration 4fs (~40 cycles). The distribution is centered around 
the line 2(w) - h - h = Si + E 2 » 5eV, with a FWHM of 
about 1.5 eV due to the finite duration of the pulse. The pulse 
(b) is a ten-cycle (~ 1 fs) sin 2 pulse. 



A. Total cross sections 

Integrating Eq. 11 over all variables including E\ and 
E2 gives, up to prefactors, the total double ionization 
cross section. The dependence on the primary photon 
energy is only implicit through the electromagnetic pulse 
entering the propagation. Within a time-dependent cal- 
culation, the resulting double ionization (DI) probabil- 
ity depends on the spectral distribution, i.e., the shape 
and duration of the laser pulse, while the fundamental 
quantity of interest, the DI cross section (DICS) at fixed 
frequency of the ionizing radiation does not. Extraction 
of the DICS therefore requires special care. 
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FIG. 2: (a) Fourier spectra of 4 fs and 1 fs sin 2 laser fields with 
a mean energy of 42 eV. (b) Double ionization probability vs. 
total energy E to t = E\ + E2 (i.e., the integral over lines with 
Etot = Ei + E2 from Fig. 1). For the long pulse, the double 
ionization probability directly reflects the Fourier spectrum. 
For the shorter pulse the electron energy distribution is in- 
fluenced by the energy-dependence of the cross section (cf. 
Fig. 8). 



For one-photon ionization, a straightforward relation- 
ship exists between the energy-dependent ionization yield 
and the cross section. From a single pulse calculation, one 
can thus obtain the cross section for all energies contained 
within the pulse [27, 48]. This is not possible without ad- 
ditional approximations for two-photon or multi-photon 
ionization, since the relation between cross section and 
yield contains an integral over intermediate energies. For 
the evaluation of this integral, the intermediate states 
and energies would have to be explicitly available. In 
the current approach, this is not easily possible without 
losing the key advantage of the time-dependent method 
of not having to construct intermediate or final states 
explicitly. 

The alternative is to use a sufficiently long pulse with 
narrow spectral width and calculate the cross section 
from the total yield with the approximation that it is 
constant over the width of the pulse. For this approxi- 
mation to be valid, the spectral width of the pulse must 
be smaller than the energy width over which the cross 
section significantly changes. We can check the conver- 
gence by varying the pulse length. Figs. 1 and 2 illus- 
trate this for both the joint two-electron energy distribu- 
tion P DC {E 1 ,E 2 ) (Fig. 1) and the integral (Fig. 2) along 
lines of constant total energy E\ + E2 in Fig. 1 for two 
different pulses: a pulse with a duration of T = 4fs con- 
taining about 40 optical cycles and one with T = 1 fs 
containing about ten optical cycles. While the 4fs pulse 
is sufficient to resolve the cross section a few eV above 
the threshold, the shorter pulse (frequently employed, see 
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refs. [25-27, 30]) results in averaging over the threshold 
region. 

Another requirement is that the pulse has to be weak 
enough such that lowest-order perturbation is applicable, 
and the ground state depletion can be neglected. We 
therefore choose a peak intensity of Iq = 10 12 W/cm 2 . 
Variation between 10 11 W/cm 2 and 10 13 W/cm 2 results 
in deviations for the total cross section at 42 eV of 
less than 0.3%. For an intensity of 10 13 W/cm 2 , the 
two-photon yield is a factor of 10 4 higher than with 
10 11 W/cm 2 . 

Another test for the applicability of perturbation the- 
ory is the linear scaling of the yield with the total dura- 
tion T of the pulse. This means that the transition rate 
must be proportional to &(t) N , where $(t) = I{t)/u is 
the photon flux and TV is the minimum number of pho- 
tons required for the process to take place. The double 
ionization yield is then given by 
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where cat is the total generalized iV-photon cross section 
for double ionization of He. Accordingly, the cross section 
is given by 
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x JJJJdE 1 dE2dn 1 dfl2P DC (E 1 ,E2,n 1 ,n 2 ), (13) 

where the effective time T e ff ,jv f° r an -ZV-photon process 
is defined as 



TeB,N — I dt 



(14) 



For a sin pulse envelope and a two-photon process, T e g 2 
is found to be 35T/128 [25, 27, 30]. Equation 13 is valid 
for direct, i.e., nonsequential, double ionization when no 
on-shell intermediate state is involved. 



B. Differential cross sections 

The triply differential cross section (TDCS) for emit- 
ting one electron with energy E\ into the solid angle Q,i, 
while the second one is emitted into Vl 2 , follows from 
Eq. 13 as 
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The TDCS presented in this paper are all calculated in 
coplanar geometry, i.e., <f>\ — <fi 2 — 0°. In the limit of an 



infinitely long laser pulse with well-defined energy (i.e., 
a delta- like spectrum), Eq. 15 becomes equivalent to 
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as calculated in time- independent approaches. Unlike the 
joint two-electron energy distribution, the TDCS as cal- 
culated by Eq. 15 is, within reasonable limits, insensitive 
to the pulse shape used in the time-dependent approach 
since the Fourier width of the pulse is accounted for by 
the integration over the energy of the second electron. 

Instead of specifying one of the energies and integrating 
over the other, it is also possible to specify energy (or 
momentum) sharing. For that purpose, we transform 
from the usual coordinates {E\,E 2 ) to (E to t,cx), with 
Etot = E\+E 2 and tan(a) = E\jE 2 . For a fixed value of 
a, the integration is performed over the total energy -Etot, 
in other words, along straight lines through the origin in 
Fig. 1. This results in the TDCS at fixed energy sharing 
(the frequently investigated case of equal energy sharing 
corresponds toa = 7r/2). 



C. Influence of final-state correlations 

Since the extraction of double ionization observables 
eventually proceeds by projection onto uncorrelated 
Coulomb final states, controlling and monitoring the ef- 
fect of residual electron-electron correlations on the cross 
section becomes important. The key point is that elec- 
tronic correlations are fully included in the initial state 
and in the time propagation and therefore in the wave 
packet at the point of projection. We monitor the resid- 
ual error by further propagating the wave function after 
the conclusion of the laser pulse (i.e., letting the electrons 
move further apart) and varying the time of projection. 
If the final state were an eigenstate of the full Hamil- 
tonian, the results would not depend on how long the 
projection is delayed. The residual dependence on the 
time of projection is thus a direct measure of the error 
introduced by the neglect of final-state correlation during 
projection. As that time is extended, this error should 
become negligible. The maximum time one can wait is 
limited in practice by the box size, as the ionized wave 
packet will hit the box boundaries at some point and be 
reflected. To test for convergence we performed one cal- 
culation with a box size of 800 a. u., using the same 4fs 
sin 2 laser pulse at 42 eV as in Fig. la, and let the wave 
function propagate for an additional 21fs after the end 
of the pulse. The doubly ionized part is still completely 
contained in the box after this time. 

Figure 3 displays the convergence of the total cross sec- 
tion as a function of the field-free propagation time r. De- 
laying the projection from r = 1 fs to r = 21 fs changes 
the total cross section by less than 2%. Extrapolating to 
infinite time (and therefore to an infinite separation of 
the two electrons, Fig. 3b) changes the cross section by 
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FIG. 3: Convergence of the total cross section with propa- 
gation time r. The cross section is calculated at different 
times t after the 4fs sin 2 pulse from Fig. 1, with intensity 
10 12 W/cm 2 (angular momentum expansion with L ma x = 3 
and 2i,max = ^2,max = 7). Plot (a) shows that while there is 
some noticeable change for short times, the results are stable 
at later times and seem to converge to a limiting value. This 
is confirmed in the inset (b), which shows the same data vs. 
1/t. Extrapolation to 1/r — > using a quadratic fit shows a 
limiting value of 0.4595, only slightly higher than the result 
obtained at r = 21 fs. Plot (c) shows the temporal evolution 
of the ratio of the expectation values of the electron-electron 
interaction energy (H12) = — f 2 1 1 ) and the total energy 
(H). 



less than 0.2%. This gives an estimate of the error due 
to projection of that order of magnitude. Furthermore, 
the electron-electron interaction energy is responsible for 
less than 1% of the total energy of the wave packet at 
t = 21 fs (Fig. 3c). The differential cross sections show 
the same qualitative convergence behavior with time as 
the total cross section. 

This shows that projecting onto products of Coulomb 
waves is not a serious limitation. In other words, when 
the electrons have had time to move apart, their in- 
teraction can be neglected when projecting onto final 
states. For higher electron energies than in this test case 
(~2.5eV), the error is expected to be even smaller and 
the convergence faster. 

Due to the fact that the coordinate space representa- 
tion of the fully correlated wave packet is available at 



the time of projection, an alternative, semi-quantitative 
check and error estimate for double ionization exists. 
From the visual inspection of snapshots of the joint radial 
distribution at different times (Fig. 4), final states rep- 
resenting double ionization can be separated from those 
representing single ionization. While the singly ionized 
part of the wave function moves parallel to the r, axes, 
the doubly ionized parts of the wave function have pos- 
itive momentum for both electrons so that they move 
away from both axes. With increasing time, the spa- 
tial overlap between singly and doubly ionized states de- 
creases, and the two contributions can be identified visu- 
ally. An estimate for an upper bound for the total double 
ionization cross section can thus be found by just inte- 
grating the radial density over the area that the doubly 
ionized wave packet occupies (7*1, r 2 > 70 a.u. in Fig. 4d). 
This integral, which still contains a small portion of single 
ionization accompanied by excitation to Rydberg states, 
gives an upper bound for the total double ionization cross 
section. In the present case (Fig. 4), the extracted es- 
timate is about 25% higher than the value determined 
by projection. This discrepancy is predominantly caused 
by the existence of high-lying Rydberg states which also 
have contributions at large values of r. 

One could expect that the choice of the final state is 
even more important when calculating differential cross 
sections because less degrees of freedom are integrated 
over. Specifically, the triply differential cross section 
(TDCS) depends on the partial-wave phase shifts, which 
may not have fully converged at the time of projection. In 
order to monitor possible errors in the angular distribu- 
tion, we have also extracted the TDCS by a complemen- 
tary method employing the coordinate representation of 
the two-electron wave packet at large propagation time, 
bypassing projection. As we are interested in the TDCS 
at equal energy sharing, we take only that part of the 
wave packet with n = ri, the part where both electrons 
have moved out to the same distance from the nucleus 
in the same time. This is what a (microscopic) time-of- 
flight detector would identify as equal-energy electrons. 
We then directly determine the angular distribution for 
this part of the wave function 
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This estimate for the TDCS, referred to in the following 
as the wave packet (WP) method, is compared with the 
projection onto Coulomb waves (Eq. 15) in Fig. 5. The 
excellent agreement we find attests to the fact that resid- 
ual errors due to final-state correlations at the point of 
projection are, indeed, negligible. Even when project- 
ing onto plane waves (not shown), we have found that 
the TDCS at equal energy sharing almost exactly agrees 
with the result obtained from projection onto Coulomb 
waves (up to a global scaling factor of about 1.1). This 
suggests that the Coulomb potential of the ionic core can 
also be neglected in the asymptotic region when only the 
differential behavior is of interest. The neglect of the 
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FIG. 4: Radial wave function density for Fig. 3, at r = 2fs (a), r = 8fs (b), r = 14 fs (c), and r = 20 fs (d) after the end of 
the pulse. For the wave function at 20 fs, the doubly ionized part of the wave function is completely contained in the box with 
fi,ra > 70 a.u. . By integrating the probability density over this region, an upper bound for the double ionization yield can 
be established. The gray lines indicate the border between singly and doubly ionized parts by visual inspection. The lighter 
gray line at ri,ra = 40 a.u. is the apparent border at r = 8fs, while the darker gray line at n,r2 = 70 a.u. is suggested by the 
distribution at r = 20 fs. The density located between the two borders contains singly ionized parts that would by mistake be 
identified as being doubly ionized at r = 8 fs due to the lack of spatial separation. 



electron-electron interaction is expected to be even less 
important. 

In addition, Fig. 5 shows the angularly resolved value 
of the electron-electron interaction energy. For this, the 
radial distance of both electrons was taken as n.,2 = 
150 a.u., which corresponds to the position of the doubly 
ionized wave packet at the point of projection. Clearly, 
the electrons only move into directions where their inter- 
action energy is negligible compared to the total energy 
of the doubly ionized wave packet. This also supports 
the finding that the electron interaction can be neglected 
when doing the projection at late times. 

The good agreement between the different methods 



used to extract total and differential cross sections can 
be understood by the following argument: The full wave 
function contains the electron correlation regardless of 
which basis it is expressed in. The only ambiguity ex- 
ists in identifying which parts of the wave packet at time 
t = t will asymptotically correspond to the situation of 
interest (i.e., double ionization in our case). If r is cho- 
sen large enough, the electrons have separated in space 
and their interaction energy is low (cf. Fig. 3c). This im- 
plies that the electrons will neither significantly deflect 
each other nor exchange energy at later times. Therefore, 
both the angular and the energy distribution are stable, 
and the momenta of the electrons at time t = r corre- 
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FIG. 5: Comparison of different methods for extracting the triply differential cross sections (TDCS) at 42 eV photon energy, 
with a 4fs sin 2 laser pulse. The data for Coulomb projection are obtained from Eq. 15 (at E\ = 2.5 eV) while the results 
labeled wave packet (WP) were obtained without transforming to momentum space (Eq. 17). The angularly resolved value 
of the electron-electron interaction energy at the position of the wave packet (n,2 ~ 150 a. u.) is also shown in comparison to 
the total energy of the doubly ionized wave packet (~ 5eV). The vertical gray line shows the ejection angle 61 of the first 
electron. The angular momentum expansion used values of L ma x = 4 and Zi.max = ta,max = 9. The radial box had an extension 
of 400 a. u., with FEDVR elements of 4a.u. and order 11. 



spond to the asymptotic momenta for t — > oo. Similarly, 
channels where one of the electrons did not gain enough 
energy to escape the Coulomb potential of the nucleus by 
the time t = r correspond to singly ionized final states, as 
the electron interaction does not provide enough energy 
to change this situation at later times. 



D. Numerical convergence tests 

In order to ensure the reliability of the calculated cross 
sections, we have performed extensive numerical testing 
and have found our results to be very well converged. 
The convergence issues addressed are (i) the radial dis- 
cretization, (ii) the temporal propagation, and (iii) the 
angular momentum expansion. 

Convergence with respect to the radial grid (i) is easy 
to achieve within the FEDVR approach. All results 
shown were obtained with finite elements of 4a.u. ex- 
tension and of order 11. Results with order 13 for el- 
ements of 4a.u. were virtually identical (within 0.02% 
for the total cross sections). Convergence of the time 
propagation (ii) using our SIL method is equally uncrit- 
ical. Even when relaxing the convergence criterion used 



for time propagation by two orders of magnitude, the 
results do not change perceptibly from those presented 
here. In addition, we also checked that our results do not 
depend on the gauge used in Eq. 2. The change in the 
total cross section when switching from velocity gauge to 
length gauge is only 0.01%. 

The final question regarding convergence concerns the 
truncation of the angular momentum expansion Eq. 6. 
As the total angular momentum L is conserved for the 
field-free Hamiltonian (because of spherical symmetry), 
the expansion does not require much higher values for 
L max than the minimum number of photons absorbed by 
the system. We have indeed found that there was no no- 
ticeable difference in any of the results between L max = 3 
or Lmax = 4. Numerically, this is especially true when 
using velocity gauge. At the low intensities used here, the 
result is well converged with L max = 3 even when employ- 
ing length gauge. The convergence with respect to single 
particle angular momenta (ii, fe), which are mixed by the 
electron-electron interaction, is much more critical. The 
size of the expansion in I2) strongly influences the ac- 
curacy of the angular distribution of the electrons and 
the degree of angular correlation. 

While the total cross section, where all angles are in- 
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FIG. 6: Convergence of the total cross section with size of 
angular momentum expansion. The total cross section is con- 
verged to within the accuracy of the method even for the 
smallest expansion in angular momenta. The differences in 
the result due to different angular basis sizes are much smaller 
than those observed when performing the projection at dif- 
ferent times after the end of the pulse (Fig. 3) (4 fs sin 2 pulse 
at 42 eV as in Fig. 1) 



tegrated over, shows almost no dependence on the size 
of the angular momentum expansion, with variations of 
less than 0.3% when (Zi, ma x, ^2,max) is increased from 
(3,3) to (9,9) (Fig. 6), a different picture emerges when 
the two-electron angular distribution is considered. The 
TDCS shows a strong dependence on the number of in- 
cluded partial waves. For the present case, convergence 
is reached when single electron angular momenta up to 
h,max = h, max = 7 are included (see Fig. 7 below). Es- 
pecially the TDCS at 9\ = 90° (where the cross section 
is very small) is very sensitive to the size of the partial 
wave expansion. 



IV. RESULTS 

In Fig. 8, we compare the present results for the total 
cross section with various published data. The calcula- 
tions were performed with a box size of 240 a. u., with 
FEDVR elements that span 4a.u. and contain 11 ba- 
sis functions. The maximum angular momentum val- 
ues are L maK — 3 for the total angular momentum and 
h,max = h, max = 7 for the individual angular momenta. 
The laser pulse envelope had a sin 2 shape, defined by 

fit) { Sin2 < < < T 

I( > ~ \ otherwise ' [LS) 

with a total duration of T — 4 fs and a peak intensity of 
Iq = 10 12 W/cm 2 . The ionization yields were extracted 
1 fs after the pulse. Following the results of section III C 
the projection error should not be larger than 2%. 

We compare our results with data from both time- 
dependent and time-independent approaches. Laulan 



and Bachau [25] solved the TDSE by means of a B- 
spline method and an explicit Runge-Kutta propaga- 
tion scheme. The double ionization probability was ob- 
tained by projecting onto uncorrelated Coulomb func- 
tions. They also included first-order correction terms in 
the representation of the double continuum (thus partly 
taking into account radial correlations). However, they 
found little difference with respect to the uncorrelated 
functions. Hu, Colgan, and Collins [26] solved the time- 
dependent close-coupling equations using finite-difference 
techniques for the spatial discretization and the real- 
space product formula as well as a leapfrog algorithm 
for temporal propagation. The double ionization prob- 
ability was also extracted by projection onto uncorre- 
lated Coulomb waves. Foumouo et al. [27] employed a 
spectral method of configuration interaction type (involv- 
ing Coulomb-Sturmian functions) and an explicit Runge- 
Kutta time propagation to solve the TDSE. The double 
continuum was generated with the J-matrix method that 
should contain angular and radial correlations to the full 
extent. In addition, they also performed calculations us- 
ing an uncorrelated representation of the two-electron 
continuum. The more recent results from Ivanov and 
Kheifets [29] are based on the time-dependent conver- 
gent close-coupling (CCC) method, taking into account 
correlations in the final state to some degree. Nikolopou- 
los and Lambropoulos [30] solved the TDSE using an 
expansion in correlated multichannel wave functions. 

Within the time-independent methods, Nikolopou- 
los and Lambropoulos [31] applied lowest-order non- 
vanishing perturbation theory (LOPT) to determine the 
generalized cross sections. Feng and van der Hart [32] 
employed i?-matrix Floquet theory in combination with 
B-splines basis sets. The data from Horner et al. [33] 
also result from LOPT calculations. They solved the 
Dalgarno-Lewis equations for two-photon absorption in 
LOPT employing exterior complex scaling (ECS) and 
also account for correlation in initial, intermediate, and 
final states. 

Overall, our results are in reasonable agreement with 
those of [25-27, 32, 33] while sizable discrepancies exist 
in comparison with those of [30, 31] as well as those of 
[27] in which corrections due to final-state correlations 
are included. Clearly, the degree of convergence of the 
present results on the few percent level as well as the up- 
per bound extracted from the radial wave packet analysis 
preclude any change of cross section by a factor of 5 — 10, 
which would be necessary to obtain values of the same 
magnitude as [27, 30, 31]. 

The experimental values of Hasegawa, Nabekawa et al. 
[34, 35] at 41.8 eV and of Sorokin et al. [36] at 42.8 eV (cf. 
Fig. 8) are compatible with most of the theoretical data. 
Due to the experimental uncertainties (e.g. the harmonic 
intensity in [34, 35] or the assumptions on the pulse shape 
and focusing conditions in [36]), the currently available 
data are not sufficient to strongly support or rule out any 
of the theoretical results. 

The present results show a more pronounced variation 
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FIG. 7: Convergence of the triply differential cross section (TDCS) with the size of the angular momentum expansion. The 
labels specify the maximum values (L maJ[ , ii, ma x, ^2,max) used in the angular momentum expansion. The vertical gray line shows 
the ejection angle #i of the first electron. The TDCS converges only for relatively large values in the angular momentum 
expansion (4fs sin 2 pulse at 42 eV as in Fig. 1). 




FIG. 8: Comparison of the total two-photon double ionization (TPDI) cross sections, obtained from Eq. 13, with T c g = 35T/128. 
In (a), the laser pulses had a sin 2 shape with a total duration of 4fs and a peak intensity of 10 12 W/cm 2 . For the results of 
Foumouo et al. [27], (NC) labels the results obtained by projecting onto uncorrelated Coulomb waves, while (FC) labels the 
results obtained using the J-matrix method, (b) shows the comparison using ten-cycle (~ 1 fs) pulses to other time-dependent 
approaches using the same pulses. Note that the results of Hu et al. [26] were rescaled by a factor of 128/70 in order to 
include the correct T c a . For both (a) and (b), the angular momenta were allowed to go up to Z/ max = 3 for the total angular 
momentum, and ?i, ma x = h,ma.x = 7 for the single electron angular momenta. The radial box had an extension of 240 a.u., with 
FEDVR elements of 4 a.u. and order 11. 
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FIG. 9: Total two-photon double ionization (TPDI) cross sec- 
tions obtained using different pulse shapes. The duration 
of the Gaussian and flat-top pulse was chosen such that the 
FWHM of the spectral distribution was identical to the one 
for the 4fs sin 2 pulse (~1.5eV). The spectral FWHM of the 
ten-cycle (~lfs) pulse is about four times larger, i.e., ~6eV. 
All other parameters were the same as for Fig. 8. 



with photon energy than other results obtained by di- 
rect integration of the time-dependent Schrodinger equa- 
tion. This can be easily explained by the fact that most 
previous work employed ten-cycle pulses. At photon en- 
ergies of 42-54 eV, this corresponds to about lfs total 
duration, and consequently, a spectral width (FWHM) 
of about 6eV (for sin 2 pulses). The results are there- 
fore an average over a rather large energy window. In 
contrast, we use pulses of 4fs duration with a narrower 
spectrum (FWHM ~1.5eV). To facilitate the compari- 
son with previous calculations we have also performed a 
calculation for a ten-cycle pulse (Fig. 8b) for which we 
find indeed better agreement. The pulse duration depen- 
dence becomes, in particular, critical near the threshold 
for sequential ionization at 54.4 eV. 

For nonsequential processes, the yield is directly pro- 
portional to the duration of the pulse, so that the cross 
section can be defined as the proportionality factor be- 
tween the double ionization rate and the photon flux $ N , 
where N is the number of photons for the direct process. 
On the other hand, the (two-photon) sequential ioniza- 
tion yield can be written as 



P™= J dt<ri*(t) J dt'a 2 <Z>(t'), (19) 

— oo t 

where o\ is the one-photon cross section for single ion- 
ization of He, and <t 2 is the one-photon cross section for 
ionization of the He + ion. Using the symmetry of the 



integrand yields 

which is proportional to the square of the total pulse du- 
ration T. Proceeding along the same lines as for Eq. 13 
by dividing the yield P^I by the pulse duration results 
in an apparent "cross section" that increases linearly with 
the pulse length, contradicting the notion of a pulse shape 
and duration independent quantity. This is not surpris- 
ing since for a two step process via on-shell intermediate 
states, a quadratic dependence on T c g (Eq. 20) is to be 
expected. If one extends the nonsequential cross section 
definition (Eq. 13) into the threshold region for the se- 
quential process, one expects a sudden rise whose height 
should be proportional to T e ff and whose width is de- 
termined by the spectral broadening of the pulse. With 
the pulses we used, the spectral width of 1.5 eV is small 
enough to observe the onset of this step discontinuity. In 
order to fully resolve the threshold behavior in a time- 
dependent calculation, even longer pulses with smaller 
bandwidth would be necessary. 

The region near the step discontinuity also provides 
a test case for the invariance of the nonsequential dou- 
ble ionization cross section under variation of the pulse 
shape. In addition to the 4fs sin 2 pulses used for most 
results shown in this paper, we also used the following 
pulse shapes: (i) a Gaussian pulse envelope and (ii) a 
flat-top pulse envelope with a sin 2 ramp on for a quarter 
of the pulse duration, constant intensity for half the pulse 
duration, and a sin 2 ramp off for the last quarter of the 
pulse. The durations of the Gaussian and flat-top pulses 
were chosen in such a way that the FWHM of the spec- 
tral distribution of the three pulse shapes was identical. 
Although all three pulses have the same spectral width, 
the distributions look different. Specifically, the spec- 
tral distribution of the flat-top pulse contains significant 
side lobes (ringing). In Fig. 9, we show that the results 
obtained for the total cross section are almost identical 
with all three pulse shapes, apart from close to the step 
discontinuity at the threshold for sequential double ion- 
ization. Note that T c g is dependent on the pulse shape 
and has to be taken into account properly. 

We turn now to the TDCS at 42 eV, the quantity 
most sensitive to the level of the underlying approxima- 
tions. The present results show qualitative agreement 
with some of the published data [26, 29], but there are 
pronounced quantitative differences. While the promi- 
nent back-to-back emission lobe (anti-)parallel to the 
laser polarization direction is well reproduced in most 
calculations (Fig. 10), the angular distribution for less 
favored emission directions (e.g. B\ = 90°) differs signifi- 
cantly from other calculations. One reason is the sensitiv- 
ity to the partial-wave expansion. In contrast to the total 
cross section, the TDCS needs a larger number of angu- 
lar momentum combinations (L, li, Z 2 ) in the expansion of 
the wave function to converge. In order to resolve angular 
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FIG. 10: Comparison of triply differential cross sections (TDCS) at 42 eV photon energy. Our data are obtained from Eq. 15 
(at Ei — 2.5 eV), with a 4fs sin 2 laser pulse. In comparison, the results of Hu et al. [26] and Ivanov and Kheifets [29] are 
shown. The vertical gray line shows the ejection angle 81 of the first electron. The angular momentum expansion used values 
of I/max = 4 and h, m ax = fe.max = 9. The radial box had an extension of 400 a.u., with FEDVR elements of 4 a.u. and order 11. 



correlations, i.e., for the triply differential cross section 
(TDCS), it is necessary to use relatively large expansions 
in single electron angular momenta. More specifically, 
good convergence of the TDCS is only reached for values 
as high as h <max = ?2,max= 7 (cf. Fig. 7), which exceeds 
the angular momentum content of most other calcula- 
tions. As discussed in section IIIC, we have alternatively 
determined the TDCS by directly analyzing the angular 
distribution of the wave packet for equal energy shar- 
ing by a radial integral constrained to equal radii. We 
find remarkably close agreement with the Coulomb pro- 
jection method (Fig. 5). The residual small deviations 
can be taken as an estimate for the uncertainty of the 
extraction method of the TDCS by Coulomb projection. 

Ivanov and Kheifets [29] take correlation in the fi- 
nal states into account using a convergent close-coupling 
(CCC) method. While the magnitude of their results is 
similar to those presented here, the shape differs consid- 
erably. In particular, they find significant probability for 
emission of both electrons in the same direction (6\ =6*2), 
where the mutual repulsion of the electrons should be 
strongest. In a very recent publication, Foumouo et al. 
[28] calculated the TDCS for equal energy sharing at 
45 eV photon energy using two different methods. The 
results obtained by projecting the final wave function on 
products of Coulomb waves resemble ours (not shown 
here for 45 eV, but the behavior is similar as for 42 eV). 



However, when correlation in the final state is taken into 
account using a J-matrix method, the results are much 
larger in magnitude (as for the total cross section, cf. 
Fig. 8) and display a shape reminiscent of the one ob- 
tained by Ivanov and Kheifets [29]. 



V. SUMMARY 

We have determined well-converged results for the to- 
tal and triply differential (generalized) cross sections for 
nonsequential two-photon double ionization of helium. 
The total cross sections agree reasonably well with a 
number of recently published papers [26, 29, 32, 33], but 
disagree with [30, 31]. While the uncorrelated results of 
[27] fit well with our data, the J-matrix results that ac- 
count for correlation, also presented in [27], are larger by 
almost an order of magnitude. In our approach, the in- 
clusion of correlation in the final double continuum states 
is bypassed by waiting long enough after the end of the 
pulse before performing the projection onto uncorrelated 
final states. 

We have also presented approximate methods to ex- 
tract both triply differential cross sections (TDCS) and 
total cross sections for nonsequential double ioniza- 
tion directly from the wave packet in coordinate space, 
thereby completely avoiding projection onto uncorrelated 
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final states. We achieve excellent agreement between 
these complementary methods providing thereby a mea- 
sure for the reliability and accuracy of the calculated 
cross sections. 

Additionally, we have analyzed the pulse length depen- 
dence of the cross sections. In most of the previous time- 
dependent approaches ten-cycle pulses have been em- 
ployed for the generalized cross sections extracted from 
the ionization yields. The resulting broad spectral width 
of the short pulse then influences the form of the cross 
section. This becomes evident from our calculations with 
considerably longer pulses. This is especially true at pho- 
ton energies above ~ 50 eV near the threshold for sequen- 
tial ionization (at 54.4 eV). Our own results for ten-cycle 
pulses, where these variations are smeared out, agree very 
well with the uncorrelated data from [27] and the results 
from [26]. We are also aware of a recent approach by 
Guan et al. [49], who obtain similar values for ten-cycle 
pulses as well. 
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